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We give a comprehensive introduction into an efficient 
numerical scheme for the minimisation of Gutzwiller en- 
ergy functionals for multi-band Hubbard models. 



Our method covers all conceivable cases of Gutzwiller 
variational wave functions and has been used success- 
fully in previous numerical studies. 
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1 Introduction In solid-state theory, multi-band Hub- 
bard models are used to study transition metals and their 
compounds. In these models only the local (atomic) part 
of the Coulomb interaction is explicitly taken into account. 
All non-local terms are included on the level of a 'Density- 
Functional Theory' calculation, which is used to set up a 
proper tight-binding Hamiltonian, see Sect. [2] 

Despite the relative simplicity of Hubbard models, 
as compared to the full electronic Hamiltonian, calculat- 
ing their properties still constitutes a very difficult many- 
particle problem. In recent years, significant progress has 
been made in this direction by the systematic study of mod- 
els in the limit of infinite spatial dimensions (D — > oo). 
The exact solution of Hubbard models in this limit leads to 
the Dynamical Mean Field Theory (DMFT), in which the 
original lattice model is mapped onto an effective single- 
impurity system that has to be solved numerically tUEEE 
[5). Although significant progress has been made in recent 
years in developing numerical techniques for the solution 
of the DMFT equations, it is still quite challenging and can 
be carried out only with limited accuracy. 

An alternative method, that also relies on infinite-Z) 
techniques, is the Gutzwiller variational approach. It al- 
lows for the approximate study of ground-state properties 
and single-particle excitations with much less numerical 
effort than within DMFT and has been applied in a number 
of works in recent years EEl[8]|9l[l0][n][I2[I3[Il[l5][l6l 
[Hl[ll?lQg20lCT A related approach 



that leads to the same energy functional for multi-band 
models is the slave-boson mean field theory [ 28 .29 . 30 31 
[321133 ! 341 . Starting from the approximate ground-state de- 
scription, it is also possible to study two-particle excita- 
tions within the 'time-dependent Gutzwiller theory' 1351 
[36][33[38][39]|40l|4Tp2l|43p4^^ . 

The main numerical problem in the Gutzwiller the- 
ory is the minimisation of the energy functional with re- 
spect to the variational parameters since their number can 
be quite large in investigations of multi-band models. We 
have developed an efficient numerical scheme for this min- 
imisation which has already been applied successfully in 
our studies on nickel ll8l fT6l and iron-pnictides [ 25 271 ■ I n 
particular, the studies on the spin-orbit coupling effects in 
nickel were numerically demanding since they required a 
rather fine energy resolution and the handling of up to 8000 
variational parameters |[T6l . To the best of our knowledge, 
no Gutzwiller minimisation of similar complexity has been 
reported in other works. We are therefore convinced that 
our minimisation algorithm will be of significant inter- 
est for all researchers who intend to apply the Gutzwiller 
theory to real materials. It is the purpose of this work to 
give detailed account of our method. Note that an alterna- 
tive method for the minimisation of a restricted class of 
Gutzwiller energy functionals has been proposed in a re- 
cent work 11501 . 

Our presentation is organised as follows. In Sections [2] 
and|3]we summarise the main results on multi-band Gutz- 
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wilier wave functions and their energy functionals in infi- 
nite spatial dimensions. Our minimisation algorithm is de- 
scribed in detail in Section [4] Some technical parts of the 
presentation are referred to four appendices. 

2 Multi-Band Hubbard models We aim to study the 
physics of multi-band Hubbard models 



H 



EE' 



er.fj' ^.f 
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(1) 



Here, we introduced the 'hopping parameters' t%'j and the 



operators c^J, which annihilate (create) an electron with 
spin-orbital index a on a lattice site i. The local Hamilto- 
nian 



H, 



E 



t;loc — c J;CTl,cr 2 l 'j )C r 1 l 'i jC r 2 



(2) 



<7"1,(72,<73,<74 



is determined by the orbital-dependent on-site energies 
£j;<Ti,cr 2 an d by the two-particle Coulomb interaction ma- 
trix elements U°~ L > (T2 < a ' 3 > a ' 4 _ We assume that the 2N spin- 
orbital states a are ordered in some arbitrary way, a = 
1, . . . , 2N where N is the number of orbitals per lattice 
site. In order to set up a proper basis of the local Hilbert 
space, we introduce the following notations for the 2 2N 
possible configurations. 

i) An atomic configuration I is characterised by the elec- 
tron occupation of the orbitals, 

7e{0;(l),...,(2J\O;(l,2),...,(2,3), (3) 
...(2iV-l,2iV);...;(l,...,2iV)}, 

where the elements in each set I = {a\,ai, . . .) are 
ordered, i.e., it is <j\ < &i < . . .. The symbol 
in (f3]l means that the site is empty. In general, we in- 
terpret the indices I as sets in the usual mathematical 
sense. For example, in the atomic configuration I\I' 
only those orbitals in / that are not in I' are occupied. 
The complement of I is I = (1,2,..., 2N)\I, i.e., in 
the atomic configuration I all orbitals but those in I are 
occupied. 

ii) The absolute value |/| of a configuration is the number 
of elements in it, i.e., 

|0| = O;|(cr 1 )| = l;|(a 1 ,<7 2 )| = 2; (4) 

...;|(l,...,2iV)|=2iV. 

iii) A state with a specific configuration / is given as 

|7> = (7+ 10> = II l°> = C • • 10) , (5) 

where the operators & a are in ascending order, i.e., it is 
<7i < era . . . < (Tin. Products of annihilation operators, 
such as 

^=n^=^-"^ ml (6) 

o-e/ 



will be placed in descending order, i.e., with <j\ > 
(T2 . • • > (Ti/i. Note that we have introduced the opera- 
tors Cj and Cj just as convenient abbreviations. They 
must not be misinterpreted as fermionic creation or an- 
nihilation operators, 
iv) The operator m/.j' = \I) describes the transfer be- 
tween configurations I' and /. It can be written as 



771/ 



r=C\C v ] I (1-rvO 



(7) 



<t"GJ 



where J = J U A special case, which derives 
from (0, is the occupation operator 



771/ 



= |I> (I\ = \{n a X{(l-n a ,) 



(8) 



The states \I) form a basis of the atomic Hilbert space. 
Therefore, we can write the eigenstates of the local Hamil- 
tonian (f2| as 

i r ) = E T ^i / ) < 9 > 

with coefficients Tj,r- With these eigenstates, the atomic 
Hamiltonian has the form 



T7j,ioc = E Ei-rrhi-r,r , 
r 



.,r,r> = in t <r'\ = Tj t rTp r , \I\ .<J'| 



(10) 



(11) 



3 Gutzwiller Energy Functional Multi-band Gutz- 
willer wave-functions have the form 



|!Pg> = Pg|^o) = n^l*°) ' 



(12) 



where tfo) is a normalised single-particle product state and 
the local Gutzwiller correlator is defined as 

Pi = J2 x i;r,r'\r) ii {r , \=J2\.p\r)u(r\ , (13) 
r,r< f 

where we introduced the matrix of variational parameters 
K-,r,r' which allows us to optimise the occupation and the 
form of the eigenstates \T)i of Pi. 

The evaluation of expectations values with respect to 
the wave function (fT2l is a difficult many -particle problem, 
which cannot be solved in general. As shown in Refs. Q 
FPU , one can derive analytical expressions for the varia- 
tional ground-state energy in the limit of infinite spatial 
dimensions (D 00). Using this energy functional for 
the study of finite-dimensional systems is usually denoted 
as the 'Gutzwiller approximation'. This approach is the 
basis of most applications of Gutzwiller wave functions 
in studies of real materials and it will also be addressed 
in this work. One should keep in mind, however, that the 
Gutzwiller approximation has its limitations and the study 
of some phenomena requires an evaluation of expectation 
values in finite dimensions 



Copyright line will be provided by the publisher 



pss header will be provided by the publisher 



3 



3.1 Local basis In general, the local density matrix 
for non-interacting electrons 

is non-diagonal with respect to <x, a'. For a fixed state \&o), 
one can always find a local basis with a diagonal density 
matrix. This will turn out to be quite useful in the minimi- 
sation with respect to the variational parameters \i-,r,r' be- 
cause, with such a basis, the energy functional has a much 
simpler form. We introduce the explicit expression of this 
simplified functional in the following Sects. 13.21 and 13.31 
If one minimises the energy with respect to |<?o}i however, 
the diagonality of (TBI is only ensured in systems with high 
symmetries. Therefore, we also need the general expres- 
sion for the variational ground-state energy with an arbi- 
trary local basis. This is given in AppendixlAl 

Note that, in general, the correlated density matrix 

is different from the non- interacting density matrix (11411 . 
In the following, however, we will frequently use the short 
term 'density matrix' for (TBI since the correlated density 
matrix ( TTBT l is not considered in this work. Moreover, we 
only study systems and wave functions which are trans- 
lationally invariant. Therefore we drop lattice site indices 
whenever this does not create ambiguities. 

3.2 Constraints As shown in Refs. f7lfT3l . it is most 
convenient for the evaluation of Gutzwiller wave functions 
in infinite dimensions to impose the following (local) con- 
straints 

(PIP)^ = 1 , (16) 

(ctptp c ff ,>* = (4 c ff ,)* . (17) 

Note that moving the operator P'P relative to cj. or c a , in 
( fTTI i does not alter the whole set of constraints. With the 
explicit form of the correlation operator (fT2l and an orbital 
basis with a diagonal local density matrix, 

C a ,a' = 5<j,a'n a , (18) 

the constraints read as 

E X r,rMr 2 < u r 2 = 1 , (19) 
r,r lt r 3 

E ^*r,r 1 ^r,r 2 m r 1 ua,r 2 ua' — &<j.<7'n<j , (20) 
r,r u r 3 
where 

\rUa)=ct\r)= £ r J>r |IUa), (21) 

m° rr , = (rhr,r'}* = E T ^ T 7,r' m ° , (22) 
i 

m^Un^Ylil-n,) . (23) 

crEl agl 

For a general orbital basis the explicit form of the con- 
straints is given in AppendixlAl 



3.3 Expectation values Each local operator 6 i; e.g., 
the local Hamiltonian (0, can be written as 

Or = E °r,r>m t ;r,r< ■ (24) 
r,r' 

In infinite dimensions, its expectation value with respect 
to (fT2l is given as 

(d) Va = £ Or 2 ,r 3 AK jA A r3>r4 m° ljr4 , (25) 
A,r 2 ,_r 3 ,A 

where the expectation values m Q r r , have been introduced 
in (l22l . Hence, the expectation value of the local Hamilto- 
nian ([Tol l becomes 

(26) 

r,r u r 3 

The expectation value for a hopping operator in infinite 
dimensions has the form 

= E & (^ 2 ) ( B h'A* 2 )#o ' (27) 

where, for an orbital basis with diagonal local density ma- 
trix, the (local) renormalisation matrix reads 

a ' A...r 4 

x((|A)(A|c ff ,)\ . (28) 

The expressions for the on-site energy and the renormali- 
sation matrix with a general orbital basis are given in Ap- 
pendixlAl 

3.4 Energy functional In a translationally invariant 
system, the expectation values, which we introduced in the 
previous section, lead to the following variational energy 
functional (per lattice site) 

E G (\ r ,r>, \%) ) = E & (qt)* E aijtT2t< ^ (29) 

(71,(72 

+ E E r^*r,r 1 ^r,r 2 m °r 1 ,r 2 ■ 
r,r x ,r 2 

Here, we introduced the tensor 

^.^^E^rCiW^ (30) 

= ^£k;, 1 ,, ! (4,„;C Mi } fl (31) 

k 

with the bare dispersion 
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The energy ([29b is a function of Xr,r' and t^o) where |i^o) 
enters d29l i, (f30b solely through the (non-interacting) den- 
sity matrix p with the elements 



Therefore, the energy 



(33) 



(34) 



has to be minimised with respect to the variational pa- 
rameters Xr,r' and the density matrix p obeying the con- 
straints (Oi (|20l), (or (O, (|77J>) and 



P 2 = P 



(35) 



This additional constraint ensures that p corresponds to a 
single-particle wave function. 

4 Numerical Minimisation of the Gutzwiller En- 
ergy Functional In principle, it is conceivable to min- 
imise the energy with respect to the variational parameters 
Xr,r' and the density matrix p simultaneously. However, 
we found it more efficient to use consecutive cycles of 'in- 
ner minimisations' (with respect to Xr,r' and with fixed p) 
and 'outer minimisations' (with respect to p and with fixed 
Xr,r') until a self-consistent minimum is reached. 

In the following we assume that all quantities in the 
energy functional and in the constraints are real. This is 
allowed since, in case of complex variational parameter or 
constraints $1% , i20\ . we may introduce the (independent) 
real and imaginary parts of these quantities. 

4.1 'Inner' Minimisation Before we explain our min- 
imisation algorithm in Sect. I4.1.2I it is essential to resolve 
the fundamental structure of our energy function. 

4.1.1 Structure of the energy function For a fixed 
density matrix p, the energy function is given as 



(Tl ,0"2 )<T-| ifTo 



Uz,z'Vzv Z ' 

Z.Z' 



(36) 



where we used the abbreviation vz for the n v variational 
parameters 

Xr ,r> 



v z = 



(37) 



which are considered as the elements of a vector v. In our 
numerical calculations we found that the inner minimisa- 
tion, as it will be described in Sect. 14.1.21 is much faster if 
we use the variational parameters d37l i instead of Xr.r'- 
The renormalisation matrix 



q° a '(v)='£SZ'(Z,Z': 



v z v z > 



(38) 



and the n c (independent) constraints (T% , (1201 1. which we 
denote as 



(Ji 



(v) = J2 M Z > Z>zvz> -gf = (1 = 1,..., nc) , 



Z.Z' 



(39) 

are quadratic functions of the variational parameters vz- 
The numbers gf in d39l correspond to the r.h.s. of Eqs. ( fT9l , 
(l20b . Note that, for a fixed density matrix p, the coefficients 
Cz,z' = {S°'{Z, Z'),fi(Z, Z'), Uz,z-} need to be calcu- 
lated only once. Moreover, we are free to work with an 
orbital basis with a diagonal local density matrix, which 
allows us to calculate these coefficients with the simplified 
energy expressions introduced in Sect. [3] It is important in 
our algorithm that the coefficients Cz,z> are stored in the 
main memory of the computer because, in this way, deriva- 
tives of all quadratic functions can be calculated very fast, 
see below. Even for large numbers n v of variational param- 
eters this can be achieved, since only a small fraction of the 
coefficients Cz.z 1 is, in fact, finite and needs to be stored. 
In case that the main-storage capacity is exceeded, there 
are several strategies to reduce the number of variational 
parameters, which we have tested. They are discussed in 
Appendix IB] 

The energy functional can be further simplified if we 
introduce the matrix 



(v) 



Z.Z' 



K (Z, Z')v z vz> 



(40) 



with the coefficients 

= K u * 2 ,a[,*> 2 SZhz,Z') . (41) 

0-2,0-2 

It allows us to write the energy as 

E G (v) - £hv)r*l(v) + Uz,z>v z v z > • (42) 



o-i,o-; 



Z.Z' 



Note that the coefficients in (1401 also need to be calculated 
only once in an inner minimisation and should be stored 
in the main memory. In this way, the energy d42l and its 
gradient E(v) with the elements 

d_ 

dv z 



E z {v) = g^E G (v) =2^2 [ E (q%(v)H&(Z,Z') 



+r%(v)S%(Z,Z')) +Uz.z'\v z 



(43) 



Z.Z' 



can be calculated very fast. The same holds for the gradi- 
ents F (v) of the constraints which have the elements 

F%(v) ^-9l(v) - 2 £ fi(Z, Z')v z > . (44) 
v z z , 

Note that in (l43l l and d44b we have used the symmetry 
Cz.z' = Cz'.z, which we are free to impose. 
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4.1.2 Algorithm for the inner minimisation We 

aim at a minimisation of the energy (l42l i in the manifold 
A4 C defined by the constraints d39b . To this end, we can 
always start our minimisation in the uncorrelated limit, i.e., 
at the point Vq (with Xr,r' = 8r,r>) f° r which t>o G M c is 
automatically fulfilled. We found numerical strategies that 
try to move exactly along M. c to be quite cumbersome. 
Therefore, starting from a certain point Vq G A4 C , we al- 
low the minimisation algorithm to violate the constraints 
by making 'short' steps to points Vi £ Ai c - To keep the 
violation of the constraints minimal, these steps have to 
take place in the subspace A4\\(vq) that is tangential to 
Ai c at the point Vq. The optimal direction of a step in 
M || (vq) is determined by the tangential component of the 
gradient E(vq) since it leads to a decrease of the energy. 
In summary, and more precisely, these ideas lead to the 
following algorithm for the inner minimisation: 

i) Find a point Vq in the variational parameter space V 
that obeys the constraints (l39l (i.e, Vq G A4 c ). 

ii) Determine the gradients F l (vo) and E(vo). 

iii) Calculate the component En (vq) of E(vq) in A4n (vq) 
by the following procedure. The gradient E(v ) is 
written as 

E(v ) = E^vo) + E ± (v ) , (45) 
where the tangential component E\\ (vq) is defined by 

E ll (v o )-F l (v o ) = VZ. (46) 

The perpendicular component can be expressed as a 
linear combination 

E ± (vo)=J^aiF l (v Q ) (47) 

i=i 

of the vectors F l (vo). In order to determine the co- 
efficients a.i, we multiply equation (|43T > with a vector 
F m (vo) and use the expansion (1471 1, This leads to 

E(v ) ■ F m (v ) = £ F l (vo) ■ F m {w) ai (48) 
i 

i 

where we used equation d46b and introduced the (sym- 
metric) matrix W(v) with the elements 

W m ,i(v) =F l (v)-F m (v) . (49) 

The linear equations d48l for ai have a unique solution, 
as long as the vectors F l (vo) are linearly independent. 
A linear dependency of these vectors can only arise if 
certain constraints d39l are redundant. In that case, the 
redundant constraints have to be eliminated right from 
the start. With the coefficients a;, we calculate the tan- 
gential component 

E ll (v ) = E(v )-J2^ Fl ( v o) ■ (5°) 

i 



of E(v ). 

iv) Make a 'proper' step in the direction of —E\\ (vo) to a 
new vector 

vi =v -!3E ll (v ) . (51) 

For the choice of the parameter (3, various strategies 
are conceivable. Since the point v\ is not in A4 C , the 
energy gain is not necessarily a useful criterion and it is 
also rather time consuming to be determined. Instead, 
we calculate 

Ag{v x ) = $>i(«i)] 2 > (52) 
i 

as a measure for the violation of the constraints and 
choose the parameter f3 such that Ag does not exceed 
a certain critical value Ag c . This critical value should 
be automatically adjusted by the algorithm to ensure 
that, after returning to the hyper-surface Ai c , there is a 
sufficient energy gain. 

v) In order to return to Ai c from the point v± £ M. c , the 
following algorithm turned out to be very useful. We 
seek a vector V\ that solves the constraint equations 
gi(vi) = and is as close as possible to v±. To this 
end, we could calculate the gradients F l (v\) and try to 
solve the set of equations 

gi(vi+Y J lmF n {v l )\=Q (53) 

^ m ' 

by a proper choice of the coefficients j m . Such an exact 
solution of equations (l53l , however, is quite time con- 
suming. Therefore, we consider the linear set of equa- 
tions 

9l(vi) + Wl,m(vihm = , (54) 
m 

which results from an expansion of (l53l to leading or- 
der in 7 m . Equations (l54l can be readily solved with 
respect to j m . This yields a new vector 

t>i v[ = vt + 7™F Tn («i) ■ (55) 

m 

which, in general, is not yet a solution of gi (v[ ) = 0. 
However, this vector is closer to M. c than V\ because 
Ag(v[) < Ag(vi). By an iteration of equations (l54l - 
$55[ we eventually approach a vector V\ G Ai c - Note 
that the fast convergence of this procedure is crucial 
for our algorithm. We have tried several other ways to 
return to A4 C that all turned out to be much slower. 

vi) If Eq(v i) < Eq{vq) we restart the procedure at point 
ii) with t>o replaced by V\. In case that Eq(vi) > 
Eq(vq), the critical value Ag c has to be lowered and 
the algorithm continues with point iv). A useful mea- 
sure for the convergence of the whole iteration is the 
norm of En. This number goes to zero near a mini- 
mum v m i a of the energy functional Eq(v) for vectors 
v G M c . 



Copyright line will be provided by the publisher 



6 



J. Bunemann et al.: Minimisation of Gutzwiller Energy Functionals 



4.2 'Outer' Minimisation With the optimum varia- 
tional parameters v mnl from the inner minimisation, de- 
scribed in Sect. 14.11 we have to minimise the energy 

Eg{p) = £ £ *Sf (p)P(j<r'),(i<r) ^ 
Z,Z' 

with respect to p. Here we introduced the renormalised 
hopping parameters 

and the renormalisation factors 

q°J(p) = Y / S:'(Z 1 Z';p)vf^r . (58) 

Z,Z' 

In addition, the (independent) constraints d76b . dTTl i, 

m(p) = £ ^ Z 'i P) v f n vz> n - 9i = (59) 

Z,Z' 

(I = 1, . . . ,n c ) , 

and (f35T > need to be obeyed. 

The local elements of the density matrix 

C<T,cr' = P(ia'),(ia) (60) 

play a special role in the energy function because only they 
enter the coefficients in (|55J, (I5B1 . ( |59l , 

Uz,z'(ft)=U z ,z'(C) , (61) 
S?(Z,Z';p) =SZ'(Z,Z';C) , (62) 
fi(Z,Z',p) =fi(Z,Z',C) . (63) 

If they are kept fixed, only the hopping term in dS6b and the 
constraint d35T l need to be taken into account in the minimi- 
sation with respect to p. This leads to a minimisation strat- 
egy which we discuss in Sect. 14.2.11 An alternative way of 
minimising (TS6b with respect to all elements of p will be 
introduced in Sect. 14.2.21 

The Hermiticity of the density matrix, ft = p, is a 
constraint which is obeyed automatically in our outer min- 
imisation algorithm in Sect. 14.2.21 To this end, however, 
the functional dependence of the energy with respect to p, 
which is not unique, must be chosen such that 

«*__(«*!_)•. m 

Op(ia), {]<?') \ a P(ja'),(ia) J 

This can always be achieved by employing the Hermiticity 
of p. We further assume that equation d64l i is also satisfied 
by the constraints d59b . 



4.2.1 Fixed local density matrix If the local density 
matrix is fixed, we have to minimise 

s G,o(p) £ E K'j P(^'U^) ( 65 > 

<*,<>' 

with respect to p obeying the constraints d35l l and ( f60b . We 
impose these constraints by means of Lagrange parameters 
i] a ^i and Q( ia .)u a i\, which leads to the 'Lagrange func- 
tional' 

Lq = Eg,o(p) ~ E ? 7 <7 > cr ' E^°'' <T ' _ P(i<r').,(ia)) 
a, a' i 

i,j cr,<r' 

As recalled in Appendix ICl the minimisation of d66b with 
respect to p leads to the effective single-particle Hamilto- 
nian 

^ ff = E + hiv«,M,„%*> ■ (67) 

The optimum single -particle state tfo) is the ground state 
of i?Q ff where the parameters r/ CT)Cr / have to be chosen such 
that C a>a > = (c], a c i a ,)^ is satisfied. 

With the state \\P) , we may determine a new ten- 
sor (f30b and start another run of the inner minimisation 
until self-consistency with respect to \\P) is reached. In 
this way, we find the ground-state energy E = Eq(C) for 
a fixed local density matrix C a , a ' . To obtain the total vari- 
ational ground-state energy, Eq (C) still needs to be min- 
imised with respect to C with the constraint of total particle 
number conservation, Co,o = N/L. Alternatively, one 
may start a self-consistency cycle of inner and outer min- 
imisation for a fixed set of 'effective crystal fields' rja;a' 
(and a fixed particle number). This defines an energy func- 
tion Eq (fj) which has to be minimised with respect to i] a .a> ■ 

Obviously, these two ways of minimising the energy 
are feasible only when the number n- s of independent ele- 
ments in C (or fields fj) is small. It can also be useful, when 
there are physical reasons to minimise Eq(C) (or Eo(fj)) 
only in some subspace of possible density matrices C (or 
fields fj)). Such a strategy has been used, e.g., in our calcu- 
lations on the spin-orbit coupling effects in nickel. There, 
we could clearly identify the relevant fields r] a : the dom- 
inant term in nickel is the effective exchange splitting ac- 
companied by a smaller orbital-energy splitting and an ef- 
fective spin-orbit coupling. In this way, the energy Eo(fj) 
had to be minimised only in a 3-dimensional subspace of 
fields fj. However, such a procedure is bound to fail when 
the number rij of parameters r\a t a' is too large and cannot 
be reduced by any physical arguments. In that case, one 
may use the algorithm which we introduce in the follow- 
ing section. 
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4.2.2 Unrestricted outer minimisation In order to 
minimise the energy with respect to all elements of the 
density matrix we impose the constraints d59l by means 
of Lagrange parameters Ai. This leads us to the functional 



U 



Eg(p) 



i 

''{}/, 



(68) 



i,j cr,cr' 



IP ~ P\(jv'),(ia) 



where Eq(p) has been defined in ( f56b . The minimisation 
with respect to p yields again an effective single-particle 
Hamiltonian of the form d67b where the fields rj a ^i are now 
given as 



dC, 



-9l(p) ■ (69) 



To determine these fields we need to calculate the Lagrange 
parameters Ai. This can by achieved if we use the fact 
that, in the variational ground state, the Lagrange func- 
tional d68l ) is also minimal with respect to the variational 
parameters Vz- This leads to the equations 



—E G (p,v) 
ov z 



which can be written in matrix- vector form as 
GA = E , 



= 
(70) 

(71) 



G 



i.z 



where G and E have the elements 
_d_ 

dv z 
d 

dv z 

The number of equations in ( FTTb is usually much larger 
then the number of parameters Ai. For physical reasons, 
however, Eq. ( |7TT i must have a unique solution. Therefore 
we can alternatively solve the equation 



-9i(p,v) 
-E G (p,v) 



(72) 



(73) 



G T GA = G T E , 



(74) 



since it gives us the same solution for A as dTIY 

Note that the calculation of the derivatives in ( |69l is 
much easier if we work with an orbital basis with a diag- 
onal density matrix, see Appendix [D] This leads us to the 
following algorithm for the outer minimisation. 

i) Set q° = Jo-.cr' and choose a reasonable set of fields 

e.g., r]^ a , ~ e<j^i with the bare on-site energies 
s (j (j' in the local Hamiltonian (fJJ. 

ii) Find the ground state of the effective Hamilto- 



If 



nian d67b with r\ a>a i = r/aja and determine C, 
Ccr j(T ' is not diagonal, find an orbital basis with a diag- 
onal local density matrix. Continue the algorithm with 
this new basis and its values for C a . a > = 5 a .a'n a and 

(J 1 ,CT 2 ,C7', ,<jL ■ 



iii) Carry out an inner minimisation, as described in sec- 
tion [4J] and determine the Lagrange parameters Ai by 
solving Eq. (l74b . 

iv)Use Eq. (|69l to determine a new set of parameters 



if^'u, ■ Set % a , = and go back to ii) until self- 



Co) 

consistency, rj a ' , 



n„„, is reached. 

la, a 



This algorithm obviously relies on a certain 'proximity' to 
the true variational ground-state, in particular, when there 
is more than one (local) minimum. In the latter case, the 
algorithm may have to be supported by a preliminary man- 
ual scan of the variational space as described in Sect. 14. 2. Tl 
Moreover, it can be necessary to introduce some kind of 
'damping' by setting 



,(0 = „(0 



(75) 



with < (3 < 1 instead of vfy^ = rj^, in step iv). The 
value of j3 must be small enough to ensure that the energy 
decreases in each step of the cycle. In our numerical tests, 
we found that f3 may sometimes have to be smaller than 
1 even in the immediate vicinity of the variational ground 
state. 

Note that the calculation of the derivatives in (|69l 
and (l72l in steps iii) and iv) of the algorithm is very much 
simplified by the fact that the local density matrix is diag- 
onal with respect to \&o). This does not mean, however, 
that the derivatives with respect to non-diagonal elements 
Ca.c' necessarily vanish, see AppendijjDl Therefore, the 
orbital basis will, in general, be changing in each cycle of 
the algorithm until a self-consistent minimum is reached. 

5 Summary In summary, we have given a detailed 
account of a numerical scheme for the minimisation of 
Gutzwiller energy functionals, which we found to be quite 
efficient in previous studies on transition metals and transi- 
tion metal compounds. We are confident that our algorithm 
is of significant interest for other researchers who intend to 
apply the multi-band Gutzwiller theory to other materials. 

A Energy functional for an arbitrary local density 
matrix The constraints (fl9l) . (l20l i for a general orbital ba- 
sis read 



^'r : r 1 ^r.r 2 ' m r 1 ,r-_ 

r,r lt r 2 



1 , 

a' Ccr,(T' j 



^*r,r 1 ^r,r 2 " l r 1 u<T,r2u 
r,r u r 2 

where 

\rua)=cl\r) = ]T Tj^lUa) 



(76) 
(77) 

(78) 



m%, = (m r ,r>)# = 2^ Tj , r T/, ^m?,,, , (79) 
IJ' 

m °i,i> = (mi.i')^a ■ (80) 
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The result for the local energy is the same as in Eq 
only with m ri r2 given by Eq. ( l79t . 

With Wick's theorem, the expectation values m° v 
in (l79l i can be written as the determinant 



(81) 



Here, Qi p are the matrices 



/ C„ 



C„ 



C' 



(82) 



in which the entries are the elements of the uncorrected 
local density matrix ( TPfl i. that belong to the configurations 



I = (tri, . . . ,<7|j|) and V 
Q J ' J in (IHTT i is defined as 



(a[, . . . , o"| 7 /|)- The matrix 



j, j 



-a 



cri,CT| 



-a 



°vi 



i - w 1J|jCT|Jj y 



wither., e J= (1,...,JV)\(/UJ'). 

The renormalisation matrix in (|27| | has the form 



3a 



E 

A,...,r 4 



(83) 



(84) 



X 



ii,U 

where the matrix Hf , contains three different contribu- 

1 1 i 1 4 

tions depending on whether the index a 1 is an element of 

h n h, h\{h n h), or J = (1, . . . , AT)\(/i U 7 4 ). With 
the abbreviation /crj = (/|c CT c CT |/) we can write Hj x l4 as 

Km =0-- U'M){h\K>\h U a')m? ij4UCT , (85) 

+ + (1 - /a',/ 4 )^?;vj 4 ) 

x(7 1 \a'|c CT ,|/ 1 ) . 

The expectation value m^^, ^ in d85l l has the same form 
as the one in (18 It . except that the index J has to be replaced 
by J\cr'. 

B Strategies to treat large numbers of 'inner' 
variational parameters Our algorithm is particularly 
fast for the inner minimisation if we can store all the 
second-order coefficients Cz,z> m the main memory of 
our computer, see Sect. 14.11 Unfortunately, this cannot al- 
ways be achieved in multi-band studies, in particular, when 
we include non-diagonal variational parameters Xrr'- hi 



this case we may try to reduce the number of variational 
parameters, e.g., by symmetry considerations, see Ap- 
pendix [bTT] Alternatively, one can employ additional nu- 
merical schemes that complement our inner minimisation 
algorithm, see Appendix lB.2l 

B.1 Reduction of the variational space It is obvi- 
ous that, due to symmetries, many parameters Xr.r' van- 
ish automatically in the variational ground state and can be 
discarded from the outset. In order to identify these param- 
eters one may use, e.g., the expectation values ( l22l which 
vanish for such parameters. 

A further reduction can be achieved if we take only 
those variational parameters into account which couple 
states \r), \r') that belong to the same (degenerate) mul- 
tiplet of the atomic Hamiltonian in ©. Such a strategy has 
been used in our calculations on the spin-orbit coupling 
effects in nickel |[T6l . Although clearly an approximation, 
this scheme is justified since one is usually bound to make 
similar approximations already on the level of the opera- 
tors in the local Hamiltonian (0. For example, in studies 
on transition metals and their compounds a spherical ap- 
proximation is often used which allows one to express all 
Coulomb-interaction parameters by the three Racah or the 
three Slater-Condon parameters. To go beyond this spheri- 
cal approximation is actually simple within the Gutzwiller 
theory, however, it increases the number of independent 
Coulomb-interaction parameters significantly. Since there 
exists no established way to calculate these parameters 
from first principles, they have to be determined by some 
fitting procedure, which only makes sense if their number 
is not too large. 

For sufficiently large Coulomb interactions, atomic 
charge fluctuations are significantly suppressed. For ex- 
ample, in elementary nickel with its approximately nine 
3d electrons per atom the occupation of states with less 
than six 3<i-electrons is negligibly small. Hence, the varia- 
tional parameters Xr.r' of such shells may be assumed to 
be diagonal or even to vanish. 

B.2 Additional numerical schemes In case that, 
even after all symmetry considerations, the number of 
variational parameters Xr r' i s st iU t0 ° large for our in- 
ner minimisation algorithm, one may employ one of the 
following numerical schemes. 

The simplest scheme is to split up the whole set of vari- 
ational parameters into sub-sets, for which the main stor- 
age of our computer is adequate and the minimisation algo- 
rithm in Sect. I4.1.2l can be applied. The minimisation with 
respect to each of these sub-sets of parameters has then to 
be repeated until a total minimum is reached. 

Another scheme is based on the observation that the 
multiplet states \r) do not necessarily have to be the eigen- 
states of our local Hamiltonian (fJJ. Instead, the states \T) 
themselves are considered as variational objects in the fol- 
lowing algorithm. 



(i) Choose a certain basis of multiplets states \T) 



(i) 
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(ii) Set \r) = \r)^ and determine the most 'relevant' 
non-diagonal variational parameters Xr,r> sucn that 
their number still allows for the use of the minimisa- 
tion algorithm in Sect. 14.1.21 A criterion for the 'rele- 
vance' of the parameters Xr,r' ma y be the size of the 
non-interacting expectation value (l22l . Alternatively 
one could use the corresponding correlated expectation 
value which can be calculated in a preceding calcula- 
tion with a diagonal variational parameter matrix A p.r- 

(iii) Determine the optimum values A^P r , of the parameters 

chosen in (ii). Calculate the eigenstates \r)(°' of the 
optimal correlation operator 



popt 



r,r> 



(86) 



(iv)Set |r)K> = |r)(°) and go back to (ii) until self- 
consistency \r)^ ~ |.r)( ) is reached. 

We have tested both numerical schemes, discussed in this 
Appendix. From these preliminary calculations, however, 
we are not yet able to draw any final conclusions on the 
efficiency of both approaches. 

C Minimisation of functions with respect to non- 
interacting density matrices We consider a general 
function E(p) of a non-interacting density matrix p with 
the elements 

Pin' = (vSk • ( 87 ) 
The fact that p is derived from a single-particle product 
wave function \<Pq) is equivalent to the matrix equation 
p 2 = p. Hence, the minimum of E{p) in the 'space' of 
all non-interacting density matrices is determined by the 
condition 

■J^-L(P) = , (88) 

"P-y'n 

where we introduced the 'Lagrange functional' 

L(p) = E(p) -J2^l,m[p 2 ~ P] m ,l 

I. Ill 

= E{p) - ^2 fy,m ( ^2 P™,pPp>1 ~ P™ l) 
l.m p 



(89) 

Pm,pPp,l ~ Pm.l I (90) 



and the matrix i? of Lagrange parameters The min- 

imisation of d89l leads to the matrix equation 



H = pfl + ttp - n 
for the 'Hamilton matrix' H with the elements 

8 



H. 



E(p). 



dp-y'n 

This equation is satisfied if p 2 = p and 

[H,p} = 0. 



(91) 



(92) 



(93) 



Hence, H and p must have the same basis of (single- 
particle) eigenvectors and, consequently, |#n) is the ground 
state of 



HI 



off 



(94) 



D Derivatives of the general energy functional 

In Sect. 14.2.21 we have to calculate the derivative of the 
ground-state energy and of the constraints with respect 
to the elements of the local density matrix, see Eq. j69l . 
Equations (l76ll-(l85ll reveal that, in fact, we only need the 
derivatives of m° v (and of rrijf- 7 ,). For a general density 
matrix C a ^i , their calculation requires an evaluation of de- 
terminants such as ( T8TT >. However, in Sect. 14.2.21 we work 
with an orbital basis for which C aM ' = 5 atCr in a . Hence 
the derivatives with respect to C at „i have a much simpler 
form. For example, for the derivatives of r , we find 

9 o x o / 1 /n<j for a el 
— — m IV = Sij-mu i (95) 
oC a .<j ' ' -1/(1 - n a ) for a f I 



for a = a 1 , and 



dC, 



- ^U'\^ {1 _ n ^( 1 _ na , ) 



(96) 



for a ^ a', where a <E I and a' 6 The derivatives of 
m^ s p are given accordingly. 
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